load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;**************************************
begin

modelname = getenv("modelname")
varname = getenv("varname")
season = getenv("season")

;modelname = "MIROC6"
;varname = "pr30d"
;season = "DJF"

Nlevel = 6                 ;6 warming levels

Nindex = 6
pct_index = (/0.90,0.95,0.98,0.99,0.995,0.999/)        ;R90/R95/R98/R99/R99.5/R99.9

nlat = 180
nlon = 360

;----------------warming periods------------------
ft := addfile("/WORK/zhangwx/EConstraint/data/CMIP6/tas/GMST/GMST_"+modelname+"_warming_levels_year.nc","r")
year_center := ft->year_center             ;(6 level) for 0,0.5,1,2,3,4C wrt baseline


;-----------------baseline threshold--------------------
diri := "/WORK/zhangwx/EConstraint/data/CMIP6/pr/"+varname+"/"+season+"/"
f0 := addfile(diri+"RX_pct_6index/RX_pct_"+varname+"_"+modelname+"_warming_levels_"+season+".nc","r")
RX_pct := f0->RX_pct(0,:,:,:)              ;(6index, lat, lon)   baseline threshold


;----------------pr data------------------
fi := addfile("/WORK/zhangwx/EConstraint/data/raw_prNd/CMIP6/"+varname+"/"+season+"/"+varname+"_"+modelname+"_historical_ssp585_1979-2100_"+season+".nc","r")
var := fi->$varname$              ;(day lat lon)   mm/day
var@_FillValue = -9999.           ;for sort purpose

time := var&time
date := cd_calendar(time,0)
year := date(:,0)


;*************************************************
exceedance_RX_pct := new((/Nlevel,Nindex,nlat,nlon/),float)
exceedance_RX_pct = exceedance_RX_pct@_FillValue

do k = 0,Nindex-1                    ;R90/R95/R98/R99/R99.5/R99.9
  do ilevel = 1,Nlevel-1

  if (ismissing(year_center(ilevel))) then
    continue
  end if

  yrstrt := year_center(ilevel)-9
  yrlast := year_center(ilevel)+10
  iselected := ind( (year .ge. yrstrt).and.(year .le. yrlast) )
  var_ilevel := var(iselected,:,:)               ;(~1800day lat lon)
  Nsample := dimsizes(iselected)

  do ilat = 0,nlat-1
    do ilon = 0,nlon-1

      idx := ind( var_ilevel(:,ilat,ilon) .ge. RX_pct(k,ilat,ilon) )
      if (all(ismissing(idx))) then
        exceedance_RX_pct(ilevel,k,ilat,ilon) = 0.
      else
        exceedance_RX_pct(ilevel,k,ilat,ilon) = tofloat(dimsizes(idx))/Nsample
      end if

    end do        ;{ilon}
  end do          ;{ilat}

  end do          ;{ilevel}
end do            ;{k}

;--------------
ilevel = 0         ;for baseline
exceedance_RX_pct(ilevel,0,:,:) = 0.10       ;R90
exceedance_RX_pct(ilevel,1,:,:) = 0.05       ;R95
exceedance_RX_pct(ilevel,2,:,:) = 0.02       ;R98
exceedance_RX_pct(ilevel,3,:,:) = 0.01       ;R99
exceedance_RX_pct(ilevel,4,:,:) = 0.005      ;R99.5
exceedance_RX_pct(ilevel,5,:,:) = 0.001      ;R99.9


;----------------------------------
exceedance_RX_pct!0 = "level"
exceedance_RX_pct&level = ispan(0,Nlevel-1,1)
exceedance_RX_pct!1 = "index"
exceedance_RX_pct&index = ispan(0,Nindex-1,1)
exceedance_RX_pct!2 = "latitude"
exceedance_RX_pct&latitude = RX_pct&latitude
exceedance_RX_pct!3 = "longitude"
exceedance_RX_pct&longitude = RX_pct&longitude

exceedance_RX_pct@long_name = "frequency exceeding baseline (1995-2014) RX_pct_threshold"
exceedance_RX_pct@units = "1"

exceedance_RX_pct&index@long_name = "90,95,98,99,99.5,99.9 percentile"
exceedance_RX_pct&level@long_name = "0C, 0.5C, 1C, 2C, 3C, 4C relative to baseline (1995-2014)"


;**********************************************************
filo := diri+"RX_pct_6index/exceedance_RX_pct_"+varname+"_"+modelname+"_warming_levels_"+season+".nc"
system("rm -f "+filo)
fo := addfile(filo,"c")
fo->exceedance_RX_pct = exceedance_RX_pct


print(varname+"_"+modelname+"_"+season+"  done")

end

